library(tidyverse)
library(glue)
library(fs)
library(raster)
library(leaflet)
# dir & csv's
#dir_dat <- "/Users/bbest/Desktop/aquamaps_ver0816c"
dir_dat <- "/Volumes/Best HD/mbon_data_big/aquamaps"
spp_csv <- file.path(dir_dat, "speciesoccursum_ver0816c.csv")
obs_csv <- file.path(dir_dat, "occurrencecells_ver0816c_head-100K-rows.csv")
#cells_csv <- file.path(dir_dat, "hcaf_species_native_ver0816c_head-100K-rows.csv")
cells_csv <- file.path(dir_dat, "hcaf_species_native_ver0816c.csv")
env_csv <- file.path(dir_dat, "hcaf_v6.csv")
spp_csv <- file.path(dir_dat, "spp.csv")
dir_spp <- file.path(dir_dat, "species_rasters")
spp <- read_csv(spp_csv) # head(occ)
#cells <- read_csv(cells_csv)
#obs <- read_csv(obs_csv)
# AquaMaps raster specifications for 0.5 degree global raster
r_na <- raster(
xmn = -180, xmx = 180, ymn = -90, ymx = 90,
resolution=0.5, crs=leaflet:::epsg4326)
r_a <- area(r_na)
r_1 <- setValues(r_na, 1)
# writeRaster(r_a, path(dir_dat, "global_area.tif"))
# writeRaster(r_1, path(dir_dat, "global_1.tif"))
env <- read_delim(env_csv, delim = "\t")
# TODO: use env.OceanArea for r_a
raster_trim <- function(r){
# source: /Users/bbest/github/nrelutils/R/raster_utils.R
# raster::trim() crazy slow compared to this
r_m <- is.na(raster::as.matrix(r))
col_notna <- which(colSums(r_m) != nrow(r))
row_notna <- which(rowSums(r_m) != ncol(r))
extent_notna <- raster::extent(
r,
row_notna[1], row_notna[length(row_notna)],
col_notna[1], col_notna[length(col_notna)])
raster::crop(r, extent_notna)
}
sp_id_to_tif <- function(sp_id){
# , cells=cells, r_na=r_na, dir_spp=dir_spp
# get cells
cells_sp <- cells %>%
filter(SpeciesID == sp_id) %>%
mutate(
cell = cellFromXY(r_na, matrix(data = c(CenterLong, CenterLat), ncol=2)))
# assign to raster
r_sp <- r_na
r_sp[cells_sp$cell] <- cells_sp$probability
r_sp <- raster_trim(r_sp) # plot(r_sp)
# write out
r_sp_tif <- glue("{dir_spp}/{sp_id}.tif")
writeRaster(r_sp, r_sp_tif, overwrite=T)
}
# calculate area of raster cells (km2)
r_a <- area(r_na)
# summarize and plot
summary(values(r_a))
plot(r_a)
cells <- read_csv(cells_csv)
spp_ids <- distinct(cells, SpeciesID) %>% pull(SpeciesID)
# takes 8.8 hours, 472 MB
for (i in 1:length(spp_ids)){ # i = 1
cat(sprintf("%04d of %d: %s -- %s\n", i, length(spp_ids), sp_id, Sys.time()))
# https://www.aquamaps.org/preMap2.php?cache=1&SpecID=Fis-22812
sp_id <- spp_ids[i]
sp_id_to_tif(sp_id)
}
# estimate time to completion ----
# 0001 of 24904: Fis-22812 -- 2018-07-01 20:49:01
# 0349 of 24904: Fis-29786 -- 2018-07-01 20:56:19
# 1214 of 24904: Fis-30623 -- 2018-07-01 21:15:14
# 24904 of 24904: WBD-Eup-82 -- 2018-07-02 05:36:41
i_all <- 24904; t_beg <- as_datetime("2018-07-01 20:49:01")
i_now <- 1214; t_now <- as_datetime("2018-07-02 05:36:41")
difftime(t_now, t_beg, units="hours") / i_now * (i_all - i_now) + t_now
# calculate avg file size ----
dir_spp <- "/Volumes/Best HD/mbon_data_big/aquamaps/species_rasters"
file_info(dir_ls(dir_spp)) %>%
summarize(
size = mean(size)) %>%
pull(size) * 24904 / 1000000
gdal2wktraster.py -r C:\Temp\TutData\SRTM\tif\*.tif -t srtm_tiled
-s 4326 -k 60x60 -I > C:\Temp\TutData\SRTM\srtm.sql
To demonstrate the methodology of “range normalization”, let’s compare a cumulative species map with and without range normalization for two species with very different ranges:
contracted: North Atlantic right whale (Eubalaena glacialis)
expansive: blue whale (Balaenoptera musculus)
For a given cell (\(c\)), we might simply take the sum of AquaMap’s Relative Environmental Suitability (RES; \(e\); ranging in value [0, 1]) per species (\(s\)) across the number of species (\(S\)) to arrive at a cumulative species richness (\(R\)) that is weighted by the species suitability within the given cell:
\[ {R_c} = \sum_{s=1}^{S} e_{s,c} \]
This approach however gives equal weight for any given cell to a species with a contracted range as an expansive one. From a population perspective, this overrepresents the expansive range. To account for small areas with high degrees of unique endemism, we can “range normalize” by dividing the given relative environmental suitability for the species and cell (\(e_{s,i}\)) by the sum of the species relative environmental suitability (\(E_s\)) across all cells (\(C\)), i.e. the whole range, and further account for cell area (\(a_c\)) and IUCN extinction risk as a species weight (\(w_s\)):
\[ E_s = \sum_{c=1}^{C} e_{s,c} * a_c \\ {RN\ R_c} = \sum_{s=1}^{S} \frac{ e_{s,c} * a_c * w_s }{ E_s } \]
Note that the (IUCN extinction risk) weight (\(w_s\)) is not part of the range normalization denominator (\(E_s\)) so when summing species those with higher extinction risk are given more weight.
# species id's for right whale and blue whale
sid_blue <- "ITS-Mam-180528"
sid_right <- "ITS-Mam-180537"
# IUCN extinction risk status
ext_blue <- "EN"
ext_right <- "CR"
# extinction risk weights
# bumping up LC to 0.2 to give credit for presence and not expecting extinct 1.2
ext_wts <- c(LC=0.2, NT=0.4, VU=0.6, EN=0.8, CR=1, EX=1.2)
w_blue <- ext_wts[ext_blue]
w_right <- ext_wts[ext_right]
get_am_raster <- function(sid, w=NULL, weight_area=T, range_normalize=T, attr_leaflet=T, extend_global=T){
# Get AquaMaps raster
# TODO: threshold param for binary output
# sid: species id # sid <- sid_right
# sid=sid_right; weight_area=T; range_normalize=T; attr_leaflet=T; extend_global=T; w=w_right
# read tif
tif <- glue("{dir_spp}/{sid}.tif")
r <- raster(tif, crs=leaflet:::epsg4326)
# set boundaries
b <- extent(r)
a = "attributes" # names(attributes(r)) # names(attributes(a))
attr(a, "bounds") <- list(lng1 = b@xmin, lat1 = b@ymin, lng2 = b@xmax, lat2 = b@ymax)
# extend to global for calculations
r <- extend(r, r_na)
attr(a, "cell_avg_e") = cellStats(r, "mean", na.rm=T)
attr(a, "sum_e") <- cellStats(r, stat='sum', na.rm=T)
if (weight_area){
r <- r * r_a
attr(a, "cell_avg_e_a") = cellStats(r, "mean", na.rm=T)
attr(a, "sum_area_km2") <- cellStats(mask(r_a, r), stat='sum', na.rm=T)
attr(a, "sum_e_area_km2") <- cellStats(r, stat='sum', na.rm=T)
}
if (range_normalize){
E_s <- ifelse(weight_area, attr(a, "sum_e_area_km2"), attr(a, "sum_e"))
r <- r / E_s
attr(a, "E_s") <- ifelse(weight_area, "sum_e_area_km2", "sum_e")
attr(a, "cell_avg_e_a_Es") = cellStats(r, "mean", na.rm=T)
}
if (!is.null(w)){
r <- r * w
attr(a, "w") <- w
attr(a, "cell_avg_w") = cellStats(r, "mean", na.rm=T)
}
attr(a, "cell_avg") = cellStats(r, "mean", na.rm=T)
if (attr_leaflet){
r_leaflet <- projectRasterForLeaflet(raster_trim(r), method="ngb")
attr(a, "r_leaflet") <- r_leaflet
}
if (!extend_global)
r <- raster_trim(r)
# append attributes
walk(names(attributes(a)), function(x){
attr(r, x) <- attributes(a)[[x]]
assign("r", r, pos = parent.frame(n=2), inherits = FALSE)})
#names(attributes(a))
#names(attributes(r))
r
}
# rasters for species relative environmental suitability [0 - 1]
r_right <- get_am_raster(sid_right, weight_area=F, range_normalize=F) # plot(r_right)
r_blue <- get_am_raster(sid_blue, weight_area=F, range_normalize=F) # plot(r_blue)
r_right_w_a_r <- get_am_raster(sid_right, w=w_right, weight_area=T, range_normalize=T) # plot(r_right)
r_blue_w_a_r <- get_am_raster(sid_blue, w=w_blue, weight_area=T, range_normalize=T) # plot(r_blue)
# names(attributes(r_right))
# attr(r_blue, "cell_avg")
# attr(r_right, "cell_avg")
# attr(r_right, "sum_area_km2")
sum_am_rasters <- function(..., attr_leaflet=T){
#add_am_rasters <- function(..., spp_wts=NULL, spp_range_normalize=T, spp_all_normalize=T){
# spp_range_normalize: for each species map, divide by sum of all pixels, before summing
# spp_all_normalize: after summing all species maps, divide by max pixel value to get range [0,1]
# spp_wts: species weights, to apply in same order as rasters. defaults to 1
# TODO: add species richness factor, so more species -> stronger result
# ...: rasters
rasters <- list(...)
s <- stack(rasters)
# if (!is.null(spp_wts)){
# stopifnot(length(rasters)==length(spp_wts))
# for (i in 1:length(spp_wts))
# s[[i]] <- raster(s, layer=i) * spp_wts[i]
# } #else {
#spp_wts = rep(1, length(rasters))
#}
# spp_wts <- ext_wts[c(ext_right, ext_blue)]
# get weighted mean
#r <- weighted.mean(s, spp_wts)
# sum across all species rasters
#r <- calc(s, function(x) )
#r <- calc(s, fun = sum, na.rm = T)
r <- sum(s, na.rm = T)
#plot(raster(s, layer=1), col=rev(rainbow(255)))
# mask
m <- overlay(s, fun=function(x) ifelse(sum(!is.na(x)) > 0, 1, NA))
r <- mask(r, m)
# trim
#r <-
#summary(r); plot(r, col=rev(rainbow(255)))
#summary(r0); plot(r0, col=rev(rainbow(255)))
#summary(r0); plot(stack, r0)
#r <- calc(s, fun = sum, na.rm = T)
#summary(r0)
# if (spp_all_normalize)
# r <- r / cellStats(r, stat='max', na.rm=T)
# project for leaflet
if (attr_leaflet)
attr(r, "r_leaflet") <- projectRasterForLeaflet(raster_trim(r), method="ngb")
r
}
#r_sum <- add_am_rasters(r_right, r_blue, spp_range_normalize=F, spp_all_normalize=F, spp_wts = NULL)
#r_nrsum <- add_am_rasters(r_right, r_blue, spp_range_normalize=T, spp_all_normalize=F, spp_wts = NULL)
r_sum <- sum_am_rasters(r_right, r_blue) # plot(r_sum)
r_nrsum <- sum_am_rasters(r_right_w_a_r, r_blue_w_a_r)
#plot(r_sum)
#plot(r_nrsum)
#r_add3 <- add_am_rasters(r_right, r_blue, spp_range_normalize=T, spp_all_normalize=T, spp_wts = NULL)
#r_add4 <- add_am_rasters(r_right, r_blue, spp_range_normalize=T, spp_all_normalize=T,
# spp_wts = ext_wts[c(ext_right, ext_blue)])
# summary(r_sum)
# summary(r_nrsum)
# summary(r_add3)
# summary(r_add4)
fit_am_raster <- function(map, r){
bounds <- attr(r, "bounds")
with(bounds, fitBounds(map, lng1, lat1, lng2, lat2))
}
add_am_raster <- function(m, r, lyr_name, r_type="spp"){
# r_type="spp" or "sum"
# r <- r_right
r <- attr(r, "r_leaflet")
max_r <- cellStats(r, stat='max', na.rm=T)
sum_r <- cellStats(r, stat='sum', na.rm=T)
rng <- c(0, 1)
if (r_type != "spp"){
rng <- c(0, max_r)
}
pal <- colorNumeric(
"inferno", rng, na.color = "transparent", reverse=T)
lab_sci <- function(v, type="numeric"){
#browser()
if (r_type == "spp"){
labs <- case_when(
v == 0 ~ "0, 0",
v == last(v) ~ as.character(glue(
"{v}, {formatC(v/sum_r, format = 'e', digits = 1)}")),
TRUE ~ "")
} else {
#if (lyr_name == "sum NR RES") browser()
labs <- case_when(
v == 0 ~ "0, 0",
v == last(v) ~ as.character(glue(
"{v}, 1")),
TRUE ~ "")
}
labs
}
m %>%
addRasterImage(
r, colors = pal, opacity = 0.6, project=F,
group=lyr_name) %>%
addLegend(
group = lyr_name, pal = pal,
position = "bottomright",
values = rng, labFormat = lab_sci,
title = ifelse(
r_type == "spp",
glue("{lyr_name} (RES, E_s)<br> sum={format(sum_r, big.mark=',', digits=1, nsmall=1)}"),
glue("{lyr_name} (max, 1)")))
}
map_am_layer <- function(lyr){
m <- leaflet() %>%
addProviderTiles(
providers$Esri.OceanBasemap, options=providerTileOptions(opacity=0.6)) %>%
add_am_raster(r_right, "right whale") %>%
add_am_raster(r_blue, "blue whale") %>%
add_am_raster(r_sum, "sum", r_type="sum") %>%
add_am_raster(r_nrsum, "sum RN", r_type="sum") %>%
addLayersControl(
baseGroups = c("ESRI Ocean Basemap"),
overlayGroups = c("right whale", "blue whale", "sum", "sum RN"),
#overlayGroups = c("right whale", "blue whale", "sum"),
options = layersControlOptions(collapsed = T)) %>%
fit_am_raster(r_blue) %>%
addMiniMap(position = "bottomleft", toggleDisplay = T, minimized = T) %>%
hideGroup("blue whale") %>%
hideGroup("right whale") %>%
hideGroup("sum") %>%
hideGroup("sum RN") %>%
showGroup(lyr)
m
}
map_am_layer("blue whale")
map_am_layer("right whale")
map_am_layer("sum")
map_am_layer("sum RN")
The equation above is problematic because it biases equatorial cells and downweights poleward cells. Another approach is to ask what does “1” mean per species richness factor?
Extinction risk (\(w_s\)): 1 = extinct, 0.8 = critically endangered, …, 0.2 = least concern
Range normalization (\(N_s\)): species normalization factor, based on the smallest species range within a given taxonomic group, eg Mammalia.
So alternatively, we describe a range normalization factor per species (\(N_s\)) that can be applied to all cells (\(C\)) regardless of the cell area (ie largest at equator, smallest at pole) relative to some minimum range size we define as 1.
\[ E_s = \sum_{c=1}^{C} e_{s,c} * a_c \\ N_s = \frac{min\{E_1, ..., E_S\}}{E_s} \\ NR_c = \sum_{s=1}^{S} e_{s,c} * N_{s} * w_s \]
Get species in Class “Mammalia” and IUCN extinction risk…
# species id's for right whale and blue whale
# sid_blue <- "ITS-Mam-180528"
# sid_right <- "ITS-Mam-180537"
library(scales)
library(rredlist)
select = dplyr::select
get_iucn_category <- function(genus_species){
# token issued to ben@ecoquants.com via http://apiv3.iucnredlist.org/api/v3/token
rl_token <- "fb963555580923eb0f0b2778060e5bcbd0cca50c0b21c416ed7ff68da9c25c11"
#genus_species <- spp_g$genus_species[1]
#genus_sp <- 'Balaenoptera brydei'
#genus_sp <- "Mesoplodon hectori"
#genus_sp <- "Mesoplodon hectori"
r <- rl_search(genus_species, key=rl_token)
# r <- rl_search_(name=genus_species, id=NULL, region=NULL, key=rl_token)
# r <- rredlist:::rr_GET(rredlist:::nir("species", "species/id", name=genus_species, id=NULL, region=NULL), id=NULL, region=NULL, key=rl_token)
# r <- rredlist:::nir("species", "species/id", name=genus_species, id=NULL, region=NULL)
#browser()
if (is_empty(r$result)){
iucn_cat <- NA
} else {
iucn_cat <- r$result$category
}
#cat(glue("{genus_species}: {iucn_cat}"), "\n")
iucn_cat
}
# bumping up LC to 0.2 to give credit for presence and not expecting extinct 1.2
iucn_weights <- c(LC=0.2, NT=0.4, VU=0.6, EN=0.8, CR=1, EX=1.2)
redo_spp_g <- F
if (!file.exists(spp_csv) | redo_spp_g){
spp_g <- spp %>%
filter(Class == "Mammalia") %>%
arrange(Genus, Species) %>%
mutate(
group = "Mammalia",
genus_species = glue("{Genus} {Species}") %>% as.character(),
iucn_category = map_chr(genus_species, get_iucn_category),
iucn_weight = iucn_weights[iucn_category],
tif = glue("{dir_spp}/{SPECIESID}.tif"),
stack_name = str_replace_all(SPECIESID, "-", "."))
write_csv(spp_g, spp_csv)
}
spp_g <- read_csv(spp_csv)
# get raster stack of species, extended to globe
s_g <- stack(sapply(spp_g$tif, function(tif) raster(tif, crs=leaflet:::epsg4326) %>% extend(r_na)))
names(s_g) <- spp_g$SPECIESID
get_spp_stat <- function(stack_name, s, stat="sum_RES"){
#stack_name <- spp_g$r[1]
r <- raster(s, stack_name)
if (stat == "sum_RES"){
cellStats(r, stat = "sum")
} else if (stat == "sum_area_km2"){
return(cellStats(mask(r_a, r), stat = "sum"))
}
}
# get stats on sum RES and area
spp_g <- spp_g %>%
mutate(
sum_RES = map_dbl(stack_name, get_spp_stat, s_g, "sum_RES"),
sum_area_km2 = map_dbl(stack_name, get_spp_stat, s_g, "sum_area_km2"),
rescale_RES = rescale(1/sum_RES),
logrescale_RES = rescale(1/log(sum_RES)))
write_csv(spp_g, spp_csv)
spp_g <- read_csv(spp_csv)
spp_g %>%
select(SPECIESID, FBname, genus_species, iucn_category, iucn_weight, sum_area_km2, sum_RES, rescale_RES, logrescale_RES) %>%
DT::datatable()
table(spp_g$iucn_category, useNA = "ifany")
##
## CR DD EN EX LC NT VU <NA>
## 2 39 13 1 47 4 10 3
iucn_weights
## LC NT VU EN CR EX
## 0.2 0.4 0.6 0.8 1.0 1.2
TODO: Range normalize across Mammalia based on latest equation using minimum species range as reference point.
library(DT)
hist(spp_g$sum_RES)
hist(spp_g$rescale_RES)
hist(spp_g$logrescale_RES)
with(spp_g, plot(sum_RES, rescale_RES))
with(spp_g, plot(sum_RES, logrescale_RES))
sum_am_stack <- function(s, spp, do_rel=T, do_ext=T, do_rng=T, do_sca=T, binary_threshold = 0.4, log_rng=T){
# parameters:
# s = # stack of species AM RES rasters for group
# spp = # table of species information
# do_rel = T # relative environmental suitability; otherwise binary
# binary_threshold = 0.4 # if do_ext=F
# do_ext = T # extinction risk
# do_rng = T # range contraction, ie endemism
# do_sca = T # final rescale [0-1]
# do_res: relative environmental suitability, otherwise binary
if (!do_rel){
v_to_binary <- function(v, na.rm=T){
if (is.na(v)) return(NA)
ifelse(v >= binary_threshold, 1, NA)
}
s <- stackApply(s, 1:nlayers(s), v_to_binary)
names(s) <- names(s)
}
# do_ext: extinction risk
if (do_ext){
s <- s * spp$iucn_weight
}
# do_rng: range contraction, ie endemism
if (do_rng){
if (log_rng){
s <- s * spp$logrescale_RES
} else {
s <- s * spp$rescale_RES
}
}
# sum
r <- sum(s, na.rm = T)
# do_sca: final rescale [0-1]
if (do_sca){
values(r) <- rescale(values(r))
}
# mask
m <- overlay(s, fun=function(x) ifelse(sum(!is.na(x)) > 0, 1, NA))
r <- mask(r, m) # plot(r);
# projectRasterForLeaflet(raster_trim(r), method="ngb") %>% mapview::mapview(attr(r, "r_leaflet"))
r
}
r_bes <- sum_am_stack(
s_g, spp_g, do_rel=F, do_ext=T, do_rng=F, do_sca=T, binary_threshold = 0.4)
r_be <- sum_am_stack(
s_g, spp_g, do_rel=F, do_ext=T, do_rng=F, do_sca=F, binary_threshold = 0.4)
r_re <- sum_am_stack(
s_g, spp_g, do_rel=T, do_ext=T, do_rng=F, do_sca=F)
r_res <- sum_am_stack(
s_g, spp_g, do_rel=T, do_ext=T, do_rng=F, do_sca=T)
r_regs <- sum_am_stack(
s_g, spp_g, do_rel=T, do_ext=T, do_rng=T, log_rng=F, do_sca=T)
r_rels <- sum_am_stack(
s_g, spp_g, do_rel=T, do_ext=T, do_rng=T, log_rng=T, do_sca=T)
map_am_sum <- function(r, name){
r_m <- projectRasterForLeaflet(raster_trim(r), method="ngb")
rng <- range(values(r), na.rm=T)
pal <- colorNumeric(
"inferno", rng, na.color = "transparent", reverse=T)
leaflet() %>%
addProviderTiles(providers$Stamen.TonerLite) %>%
addRasterImage(
r_m, colors = pal, opacity = 0.8, project = F, group = name) %>%
addLegend(
group = name, pal = pal,
position = "bottomright",
values = rng,
title = name)
}
map_am_sum(r_bes, "Mammalia: Binary<br> * Extinction")
map_am_sum(r_res, "Mammalia: Relative<br> * Extinction")
map_am_sum(r_regs, "Mammalia: Relative<br> * Extinction<br> * Endemism")
map_am_sum(r_rels, "Mammalia: Relative<br> * Extinction<br> * log(Endemism)")
Compare species ranges with available OBIS data to determine degree of coverage, especially for determining where we expect concentrations of endangered species and yet have no observational data.
For example with North Atlantic Right Whale (Eubalaena glacialis):
OBIS observations: 5695 records for Eubalaena glacialis
The topic of species home ranges is rich. Am keeping track of references in the mbon-indicators Zotero group, eg folder species home ranges:
Peterson, A. T. (2017). Problems with reductive, polygon-based methods for estimating species’ ranges: reply to Pimm et al. 2017: Estimating Species’ Ranges. Conservation Biology, 31(4), 948–951. https://doi.org/10.1111/cobi.12929
Cochrane et al (2016) What Is Marine Biodiversity? Towards Common Concepts and Their Implications for Assessing Biodiversity Status Frontiers in Mar Sci
Pimm et al (2014) The biodiversity of species and their rates of extinction, distribution, and protection Science
rpostgis::: R Interface to a ‘PostGIS’ Database
pgWriteRast(conn, name, raster, bit.depth = NULL, blocks = NULL, constraints = T, overwrite = F): write raster to PostGIS database table
pgGetRast(conn, name, rast = "rast", bands = 1, boundary = NULL): load raster from PostGIS databasepgListRast(conn): list geometries/rasterspgGetBoundary(conn, name, geom = "geom", clauses = NULL): retrieve bounding envelope of geometries or rastersst_intersection() function and st_intersects() operator which operates directly on raster and geometries by vectorizing only the necessary part of the raster before doing the intersection and is one of the main feature of PostGIS WKT Raster”docker run -d -p 80:80 -h cartodb.localhost sverhoeven/cartodbUsing Mapbox Vector Tiles in CARTO for Maps & Location Apps — CARTO Blog
raster2pgsql
ST_Shift_Longitude - Exploring Data Formats and Fields
Here’s how the presence of species (binary, from IUCN range map or AquaMaps RES > 0.4) in a given region was assessed for average extinction risk (using IUCN extinction category.)
\[ x_{SPP} = \frac{\sum_{k=1}^{M} (\frac{\sum_{i=1}^{N} w_{i}}{N}) * A_{C}}{A_{T}} \]
\(M\) = number of grid cells in the assessment region
\(N\) = number species in a grid cell \(c\)
\(A\) = total area of a grid cell [\(c\)] the assessment region [T]
\(w_{i}\) = status weight assigned per threat
assessed species list and maps: IUCN